NASA / TM-2000-210552 



Exploratory Studies in Generalized 
Predictive Control for Active Aeroelastic 
Control of Tiltrotor Aircraft 


Raymond G. Kvaternik and Jer-Nan Juang 
Langley Research Center, Hampton, Virginia 

Richard L. Bennett 

Bell Helicopter Textron, Inc., Fort Worth, Texas 


October 2000 



The NASA STI Program Office ... in Profile 


Since its founding, NASA has been dedicated 
to the advancement of aeronautics and space 
science. The NASA Scientific and Technical 
Information (STI) Program Office plays a key 
part in helping NASA maintain this important 
role. 

The NASA STI Program Office is operated by 
Langley Research Center, the lead center for 
NASA's scientific and technical information. 
The NASA STI Program Office provides 
access to the NASA STI Database, the largest 
collection of aeronautical and space science 
STI in the world. The Program Office is also 
NASA's institutional mechanism for 
disseminating the results of its research and 
development activities. These results are 
published by NASA in the NASA STI Report 
Series, which includes the following report 
types: 

• TECHNICAL PUBLICATION. Reports 
of completed research or a major 
significant phase of research that 
present the results of NASA programs 
and include extensive data or theoretical 
analysis. Includes compilations of 
significant scientific and technical data 
and information deemed to be of 
continuing reference value. NASA 
counterpart of peer-reviewed formal 
professional papers, but having less 
stringent limitations on manuscript 
length and extent of graphic 
presentations. 

• TECHNICAL MEMORANDUM. 
Scientific and technical findings that are 
preliminary or of specialized interest, 
e.g., quick release reports, working 
papers, and bibliographies that contain 
minimal annotation. Does not contain 
extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. 
Collected papers from scientific and 
technical conferences, symposia, 
seminars, or other meetings sponsored 
or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, 
often concerned with subjects having 
substantial public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign 
scientific and technical material 
pertinent to NASA's mission. 

Specialized services that complement the 
STI Program Office's diverse offerings 
include creating custom thesauri, building 
customized databases, organizing and 
publishing research results ... even 
providing videos. 

For more information about the NASA STI 
Program Office, see the following: 

• Access the NASA STI Program Home 
Page at http://www.sti.nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA STI 
Help Desk at (301) 621-0134 

• Phone the NASA STI Help Desk at 
(301) 621-0390 

• Write to: 

NASA STI Help Desk 

NASA Center for AeroSpace Information 

7121 Standard Drive 

Hanover, MD 21076-1320 



NASA / TM-2000-210552 



Exploratory Studies in Generalized 
Predictive Control for Active Aeroelastic 
Control of Tiltrotor Aircraft 


Raymond G. Kvaternik and Jer-Nan Juang 
Langley Research Center, Hampton, Virginia 

Richard L. Bennett 

Bell Helicopter Textron, Inc., Fort Worth, Texas 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virginia 23681-2199 


October 2000 



This paper was presented at the American Helicopter Society Northeast Region Specialists' Meeting on 
Active Controls Technology, Bridgeport, CT, October 4-5, 2000. The full paper, which consisted of a set 
of charts and a text explanation of each chart, is reproduced in this report. 


The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 
official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics 
and Space Administration. 


Available from: 

NASA Center for AeroSpace Information (CASI) 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301)621-0390 


National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, VA 22161-2171 
(703) 605-6000 






Exploratory Studies in Generalized Predictive Control 
for Active Aeroelastic Control of Tiltrotor Aircraft 


Raymond G. Kvatemik and Jer-Nan Juang 
NASA Langley Research Center 
Hampton, VA 

and 

Richard L. Bennett 
Bell Helicopter Textron, Inc 
Fort Worth, TX 


Abstract 

The Aeroelasticity Branch (AB) at NASA Langley Research Center has a 
long and substantive history of tiltrotor aeroelastic research. That research has 
included a broad range of experimental investigations in the Langley Transonic 
Dynamics Tunnel (TDT) using a variety of scale models and the development of 
essential analyses. Since 1994, the tiltrotor research program has been using a 
1 1 5-scale, semispan aeroelastic model of the V-22 designed and built by Bell 
Helicopter Textron Inc (BHTI) in 1983. That model has been refurbished to form 
a tiltrotor research testbed called the Wing and Rotor Aeroelastic Test System 
(WRATS) for use in the TDT. In collaboration with BHTI, studies under the 
current tiltrotor research program are focused on aeroelastic technology areas 
having the potential for enhancing the commercial and military viability of 
tiltrotor aircraft. Among the areas being addressed, considerable emphasis is 
being directed to the evaluation of modern adaptive multi-input multi-output 
(MIMO) control techniques for active stability augmentation and vibration 
control of tiltrotor aircraft. As part of this investigation, a predictive control 
technique known as Generalized Predictive Control (GPC) is being studied to 
assess its potential for actively controlling the swashplate of tiltrotor aircraft to 
enhance aeroelastic stability in both helicopter and airplane modes of flight. This 
paper summarizes the exploratory numerical and experimental studies that were 
conducted as part of that investigation. 
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Some Relevant Background 


AB/TDT has a long history of tiltrotor aeroelastic research 
Research has included a broad range of experimental 
investigations and the development of essential analyses 
Since 1994, tiltrotor research program has been using a 
1/5-scale, semispan aeroelastic model of the V-22 

AB/BHTI cooperative study of promising tiltrotor 
aeroelastic technology areas underway since 1994 

Evaluation of modem adaptive MIMO control techniques 
for active stability augmentation a major activity 

GPC-based method for actively controlling swashplate to 
enhance stability under investigation since 1997 


The Aeroelasticity Branch (AB) at NASA Langley Research Center has 
a long and substantive history of tiltrotor aeroelastic research (ref. 1). That 
research has included a broad range of experimental investigations in the 
Langley Transonic Dynamics Tunnel (TDT) using a variety of scale models 
and the development of essential analyses. Since 1994, the tiltrotor research 
program has been using a 1/5-scale, semispan aeroelastic model of the V-22 
designed and built by Bell Helicopter Textron Inc (BHTI) in 1983. That model 
has been refurbished to form a tiltrotor research testbed called the Wing and 
Rotor Aeroelastic Test System (WRATS) for use in the TDT. hi collaboration 
with BHTI, studies under the current tiltrotor research program are focused on 
aeroelastic technology areas having the potential for enhancing the commercial 
and military viability of tiltrotor aircraft. Among the areas being addressed, 
considerable emphasis is being directed to the evaluation of modem adaptive 
multi-input multi-output (MIMO) control techniques for active stability 
augmentation and vibration control of tiltrotor aircraft. As part of this 
investigation, a predictive control technique known as Generalized Predictive 
Control (GPC) is being studied to assess its potential for actively controlling 
the swashplate of tiltrotor aircraft to enhance aeroelastic stability in both 
helicopter and airplane modes of flight. This paper summarizes the exploratory 
numerical and experimental studies that were conducted as paid of that 
investigation. 
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Approaches to Control 


Conventional control approach 

— System modeling and verification (off-line) 

— Controller design and verification (off-line) 

Adaptive control approach 

— On-line/real-time system identification 

— On-line/real-time controller design 

— User specification of model order, horizons, and control weights 

Autonomous control approach 

— System identification/controller design combined into one step 

— Automated estimation of system order and disturbances 

— Automated adjustment of control weights and control gains 

— All done on-line and in real-time 


The approaches to active control may be grouped into three general 
categories: conventional control, adaptive control, and autonomous control. 

Conventional approach: Uses an explicit plant model represented by a 
transfer function or state-space model. Determination and verification of the 
system’s transfer function and the design of the controller are done off-line. 
LQR, Hoc, H2, and p-synthesis are examples of such approaches. 

Adaptive approach: The term adaptive as used here refers to a real-time 
digital system that updates the system identification (SID) parameters based on 
input/output data, uses the updated system parameters to compute a new set of 
control parameters, and then computes the next set of commands to be sent to 
the actuators, with all computations being done on-line. The user specifies the 
order of the system model, the prediction and control horizons, and the control 
weights. 

Autonomous approach: SID and controller design are combined into one 
step. Estimation of system order and disturbances, and adjustment of control 
weights and control gains (e.g., via fuzzy logic, neural nets, and genetic 
algorithms) are all automated. All computations are done on-line and in real- 
time. Autonomous control is a very active area of research at this time. 

The work addressed in this paper falls into the second category. 
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Generalized Predictive Control 


Belongs to class of model-based predictive controllers used 
in process control industries since late 1970s 
Method introduced in 1987 by Clarke of Oxford University 
Input-output map of system used to form multi-step output 
prediction equation over a finite prediction horizon while 
subject to controls imposed over a finite control horizon 
Optimal control sequence determined by minimizing 
deviation of future system outputs from desired (target) 
response, subject to penalty on control power 

A novel version of GPC algorithm developed at NASA- 
Langley in 1997 for efficient computation and unknown 
disturbance rejection 


Predictive controllers introduced in the chemical industries for controlling 
chemical processes have found applications in a wide variety of industrial processes 
(e.g., ref. 2). Predictive control refers to a strategy wherein the decision for the 
current control action is based on minimization of a quadratic objective function that 
involves a prediction of the system response some number of time steps into the 
future. A variety of predictive controllers have been proposed (e.g., ref. 3). Among 
these, Generalized Predictive Control (GPC), which was introduced in 1987 (refs. 4- 
5), has received notable attention by researchers. GPC is a time-domain multi- 
input-multi-output (MIMO) predictive control method that uses a finite-difference 
representation for the input-output map of the system. The input-output equation is 
used to form a multi-step output prediction equation over a finite prediction horizon 
while subject to controls imposed over a finite control horizon. The control to be 
imposed at the next time step is determined by minimizing the deviation of the 
predicted controlled plant outputs from the desired (or target) outputs, subject to a 
penalty on control effort. 

A novel version of the GPC procedure was developed at NASA Langley 
Research Center in 1997 for efficient computation and unknown disturbance 
rejection by Dr. Jer-Nan Juang and his associates. Their work has resulted in a suite 
of MATLAB m-files that have been collected into a Predictive toolbox that can be 
used by researchers for GPC studies. A summary of the SID and control theory 
underlying their development is found in references 6-13, among others. These 
references were the primary sources of information for the work reported in this 
paper. 


4 




Adaptive Control Process 


Disturbance (d) 



The essential features of the adaptive control process used in the present 
GPC investigation are depicted in the diagram above. The system (plant) has r 
control inputs u, m measured outputs y, and is subject to unknown external 
disturbances d. Measurement noise is also present. There are two fundamental 
steps involved: (1) identification of the system; and (2) use of the identified 
model to design a controller. A finite-difference model in the form of an Auto- 
Regressive moving average with exogenous input (ARX) model is used here. 
This model is used for both system ID and controller design. System 
identification is done on-line in the presence of any disturbances acting on the 
system, as indicated in the center box of the diagram, hi this way, an estimate of 
the disturbance model is reflected in the identified system model and does not 
have to be modeled separately. This approach represents a case of feedback with 
embedded feedforward. Because the disturbance information is embedded in the 
feedforward control parameters, there is no need for measurement of the 
disturbance signal (ref. 13). The parameters of the identified model are used to 
compute the predictive control law. A random excitation u id (sometimes called 
dither) is applied initially with u c equal to zero to identify the open-loop system. 
Dither is added to the closed-loop control input u c if it is necessary to re-identify 
the system while operating in the closed-loop mode. 
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Form of Model and Control Law Equations 

ARX model used for system identification 
y(k) = aj(k-l) + a 2 y{k-l) + • • • + a p y(k - p) 

+ p n u(k) +J3 l u(k- 1 ) + ••• +J3u{k-p ) 

Control law equation 

u c {k) = a[y(k-l) + a c 2 y(k- 2 ) + • • • + a c p y(k - p) 
+fiu(k- 1 ) +fcy{k- 2 ) + ... +F p u{k-P ) 


The relationship between the input and output time histories of a MIMO 
system are described by the time-domain AutoRegressive exogenous (ARX) 
finite-difference model shown in the first equation. This equation states that the 
current output y(k) at time step k may be estimated by using p sets of the 
previous output and input measurements, y(k-l y(k-p) and u(k-l u(k-p); 
and the current input measurement u(k). The integer/? is called the order of tire 
ARX model. The coefficient matrices a, and />, appearing in this equation are 
referred to as observer Markov parameters (OMP) or ARX parameters and are 
the quantities to be determined by the identification algorithm. Closed-loop 
robustness is enhanced by performing the system identification in the presence 
of the external disturbances acting on the system, thereby ensuring that 
disturbance information will be incorporated into the system model. The goal of 
SID is to determine the OMP based on input and output data. The OMP may be 
determined by any SID techniques that returns an ARX model of the system. 

The ARX model is used to design the controller and leads to a control 
law that in the case of a regulator problem has the general form given by the 
second equation. This equation indicates that the current control input u c (k) may 
be computed using p sets of the previous input and output measurements. The 
coefficient matrices a j and fi appealing in this equation are the control gam 
matrices. The derivation of the parameters appearing in the SID and control law 
equations is described below. 
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System Identification 

• The digitized input and output time histories at / time 
points are used to form the data matrices y and V in 
the equation — 

y = YV 

according to 

y = [y(0) y(l) y(2) ■■■y(p) •••y(/-l)] 

rn X / 


u(0) i/(l) 

u(2) ■ 

• u(p) 

//(/ — 1) 

v(0) 

v(l) • 

■ v(p-l) 

- v(/-2) 


v(0) • 

• v(p- 2) 

- v(/-3) 



v(0) 

••• v(l-p- 1) 


[/• + (/■ + m)p]xl 


System identification in the presence of the operational disturbances 
acting on the system is the fust of the two major computational steps. The 
external disturbances acting on the system are assumed to be unknown 
(unmeasurable). The number of control inputs is r and the number of measured 
outputs is m. The order of the ARX model (p) and the number of time steps (/) 
must be specified. Some guidelines for their selection are given later. 

The system is excited with band-limited white noise for SID. These 
random excitations are input to all r control inputs simultaneously and the m 
responses measured. The digitized input and output time histories (u and y) at / 
time points are then used to form the data matrices y and V as indicated in the 
slide. The sizes of the resulting arrays are noted. The equation shown follows 
from writing the discrete-time state-space equations for a LTI system at a 
sequence of time steps k = 0, 1, ..., I - 1 and grouping them into the above matrix 
form. The vector v is defined on the next slide. 

hi forming the matrices above, it has been assumed that the state matrix A 
is asymptotically stable so that for some sufficiently large p, A i: ~ 0 for all time 
steps k > p, and that a observer has been added to the system. It is through these 
expedients that the matrix V is reduced to a size amenable for practical 
numerical computation of its pseudo-inverse. The SID process yields OMP 
rather than system Markov parameters (SMP) because of the inclusion of an 
observer. A complete discussion of these aspects of the development may be 
found in reference 6. 
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System Identification (Concluded) 


The vector v(k) in the matrix V is formed according to 

(r + m)X 1 

Y is the matrix of observer Markov parameters (OMP) 
which are to be identified and has the form 


: A, A 


a , 








vj ^3 § p 

mxr mxr mxm mxr mxm mXr mXm mXr mXm 
mx[r + (r + m)p] 

The solution for ^ is given by 

F= y v f =y v r [v V 7 ] 1 

where f denotes the pseudo-inverse. 


The vector v(k) that appears in the data matrix V is formed from the 
vectors u(k) and y(k) as indicated above. Because of the introduction of an 
observer into the discrete state equations for the system (ref. 6), the quantities a, 
and Pi in Y are referred to as observer Markov parameters (OMP). The 
solution for Y , the vector of observer Markov parameters, is obtained by solving 
the equation y = Y * V for Y as shown in the bottom equation in the slide 
above. If the product W l is a well-conditioned matrix, the ordinary inverse can 
be taken as shown. Otherwise, a pseudo-inverse must be used. It should be 
noted that because the size of VV l is much smaller than V, a pseudo-inverse may 
be appropriate even if the product is well conditioned. 
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Derivation of Multi-Step Output Prediction Equation 

The input-output map is given by the ARX equation 
y(k) = a l y(k- 1) + a 2 y(k- 2) + ••• + a p y(k-p) 

+ j3 0 u(k) + fiyuik- 1) + ••• +J3 p u(k-p) 

Using this equation, the output at time step k+j can be written as 
y(k + j) = a[ i) y(k~\) + a ( 2 ]) y{k-l) + ••• +a { p j) y(k-p) 

+ fio u (k + j) + P<pu(k + j - 1 ) + ••• + Po^uik) 

+P[ ,) u(k-X)+P2 ) u(k - 2 ) + ••• +/3pJ\i(k-p) 

where the coefficient matrices are combinations of a x and 


The one-step ahead ARX equation describing the input-output map of the 
system is the starting point for deriving the multi-step output prediction equation 
that is needed for designing the MIMO controller. Using the one- step ahead 
equation shown at the top of the slide, the output at time step k+j may be written 
in the form shown in the second equation, where the coefficient matrices are 
given by recursive expressions involving the matrices a, and /i, appearing in the 
ARX equation (ref. 10). The matrices a, and />, are determined by the system 
identification process described earlier. 

The second equation states that the output y(k+j) at time step k+j may be 
estimated by using p sets of the previous output and input measurements, 
y(k-l y(k-p) and u(k-\),...,u(k-p), and the (unknown) current and future 
inputs, u(k), u(k+\),...,u(k+j). Letting / range over the set of values j = 1, 2, 
h p - 1, the resulting equations can be assembled into a multi-step ahead output 
prediction equation having the form indicated in the next slide. 
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Derivation of Prediction Equation (Concluded) 

Letting j range over the set of values j = 1 , 2, . . . , q, . . . , h - 1 , the resulting 
equations can be assembled into the multi-step output prediction equation: 


Current & 

_ Current & 

future Output 

| f uture Input 

y K (k) = , T u h (k)+ 

h„r x 1 


B 


Past Input 

I 


U, 


Past Output 

I 


h p m x 1 


h p m x h c r 


h m x pr 


( k-p)+ A yJk-p ) 


pr x 1 


h m x pm 


pm xl 


where the coefficient matrices 'T , B, A are formed from combinations 
of the observer Markov parameters Ci t and f3 t appearing in the ARX 
equation describing the input-output map of the system, and where h p is 
the prediction horizon, h c is the control horizon, and h c C h p . 


GPC is based on system output predictions over a finite horizon h p known 
as the prediction horizon. To predict future plant outputs, some assumption 
needs to be made about future control inputs. In determining the future control 
inputs for GPC, it is assumed that the control is applied over a finite horizon h c 
known as the control horizon. Beyond the control horizon, the control input is 
assumed to be zero. In GPC, the control horizon is always equal to or less than 
the prediction horizon. 

Tire general form of the multi-step ahead prediction equation indicating the 
sizes of the constituent matrices is shown in this slide. The coefficient matrices 
T, B, and A are formed from combinations of the observer Markov parameters 
a, and Pi appearing in the ARX equation describing the input-output map of the 
system. The objective is to predict the output for h p time steps ahead, given the 
input for h c steps ahead, where h c ^ h p . The quantity yh p (k) is the vector 
containing the predicted future outputs, whereas UhJk) is the vector containing 
the future control inputs yet to be determined. The quantities u v (k-p) and y p (k-p) 
are vectors containing the previous p sets of control inputs and outputs, 
respectively. 
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Expanded Form of Multi-Step Output Prediction Equation 
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The expanded form of the multi-step output prediction equation for the 
case in which the control horizon is less than the prediction horizon is shown in 
this slide. The OMP a, and fi, determined in system identification form the first 
block rows of the coefficient matrices 7 ~, A. , and 3 in the equation. The terms 
in the remaining rows are computed using the recursive relations indicated in the 
boxes (ref. 10). All terms in the above equation are known, except for the h c 
sets of future commands, and the h p sets of predicted responses. The goal of the 
GPC control algorithm is to determine the set of future commands u(k), u(k+ 1), 
..., u(k+h c - 1 ) that are required to achieve a desired predicted response y(k), 
y{k+\), ...,y(k+h p -\). 

It should be remarked that the system Markov parameters (SMP), which 
are commonly used as the basis for identifying discrete-time models for linear 
dynamical systems, form the first block column in the matrix l ; the remaining 
block columns are formed from subsets of the SMP. The Markov parameters 
are the pulse response of a system and are unique for a given system. The 
discrete-time state- space matrices A, B, C, and D are embedded in the SMP. 

The multi-step output prediction equation is used to define an objective 
function whose minimization with respect to Uhc(k) leads to the control law from 
which a vector of future control inputs can be computed using the p sets of 
previous control inputs and measured outputs. 
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Derivation of Control Law 


Define an error function s: 

£ = y T ( k ) - y hp (k ) = yr (*) - T u h c (*) - Bu P ( k ~p)~ Ay p (k - p) 

Define an objective function J: 

J =£ T R£ + u h T Qu hc 

Minimize J with respect to u h (k ) : 

u hc (k) = ~ ( rT RT + Q) ] tT r (-Jr ( k) + Bu p (k-p)+Ay p (k-p ) ) 

h c rx 1 

Retain the first component (first r rows) of u i, r (k) : 
u c (*) = - fyA k ) + ( k ~ p)+ aC y„ ( k - p) 

rx 1 


The predictive control law is obtained by minimizing the deviation of the 
predicted controlled response (as computed from the multi-step output 
prediction equation) from a specified target response over a prediction horizon 
h p . To this end, one first defines an error function that is the difference between 
the desired (target) response yr (k) and the predicted response }'hp(k) and forms 
an objective function J quadratic in the error and the unknown future controls. 
Two weighting matrices are included in the objective function: Q (symmetric 
and positive-definite ) is used to limit the control effort and stabilize the closed- 
loop system; R (symmetric and positive- semidefinite) is used to weight the 
relative importance of the differences between the target and predicted 
responses. Typically, Q and R are assumed to be diagonal and for Q to have the 
same value w c along its diagonal and R to have the same value viy along its 
diagonal. Minimizing J with respect to Uhc(k) and solving for Uh c (k) gives the 
control sequence to be applied to the system over the next h c time steps. The 
first r values (corresponding to the first future time step) are applied to the r 
control inputs, the remainder are discarded, and a new control sequence is 
calculated at the next time step. 

The target response is zero for a regulator problem and non-zero for a 
tracking problem. Q must be tuned to ensure a stable closed-loop system. 
Typically, h c is chosen equal to h p . However, h c may be chosen less than h p 
resulting in a more stable and sluggish regulator. 
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Features Characterizing Method 


Predictive controllers employ a strategy wherein current 
control actions are based on prediction of system response 
at a number of time steps into the future 
Approach allows one to work directly with an input-output 
(ARX) model of the system 

SID process employs an observer to enable numerical 
computation of pseudo-inverse for calculation of OMP 

Controller is thus inherently observer-based but no explicit 
consideration of its presence needs to be taken into account 

System identification is performed in the presence of any 
disturbances acting on system 


Predictive controllers such as GPC employ a strategy wherein the values of 
the current control inputs are based on the predicted responses at a number of 
time steps into the future. Such an approach lends itself to the use of an ARX 
model to represent the system for both system identification and controller 
design. The SID process used makes recourse to an observer to enable 
numerical computation of the pseudo-inverse needed for calculation of the OMP 
that comprise the coefficients of the ARX equation. The controller is thus 
inherently observer-based but no explicit consideration of the observer needs to 
be taken into account in the implementation. 

hi practice, the disturbances acting on the system are unknown or 
unmeasurable. However, as discussed in reference 13, by performing the SID in 
the presence of the external disturbances acting on a system, a disturbance 
model is implicitly incorporated into the identified observer Markov parameters. 
However, the identified model must be larger than the true system model to 
accommodate the unknown disturbances. 
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Features Characterizing Method (Continued) 


Disturbance model is implicitly incorporated into the 
identified coefficients (OMP) of the ARX model 


Effects of unknown disturbances embedded in coefficient 
matrices A and B (and hence« c and P c ) because SID 
done with disturbances acting on system 


Estimated order of ARX model used for system ID given 
by 


p > ceil 


rt + 11 

system slates disturbance states 

m 




# Prediction and control horizons are set according to 
h p >p h c < h p 


The disturbance model is implicitly incorporated into the identified 
coefficients (OMP) of the ARX model. Thus, the effects of the unknown 
disturbances acting on the system are embedded in the matrices A and B , and 
hence the control law matrices of and [> . 

An expression for estimating the order of the ARX model that is to be used 
for SID is given in the slide. The number of system states is typically chosen to 
be twice the number of significant structural modes; the number of disturbance 
states is set to twice the number of frequencies in the disturbance; m is the 
number of output measurements. If measurement noise is of concern, the order 
of p so computed should be increased to allow for computational poles and zeros 
to improve system identification in the presence of noise, hi practice, simply 
choosing p to be 5-6 times the number of significant modes in the system is 
often adequate. 

The prediction and control horizons are set according to the relations 
indicated at the bottom of the slide. Although h p can be set equal to p, h p is 
typically set greater than p to help limit control effort so as not to saturate the 
control actuators. If the control horizon is greater than the system order a 
minimum energy (minimum norm) solution is obtained wherein the commanded 
output is shared so that the control actuators don’t fight each other. If h p is set 
equal to p one obtains a so-called deadbeat controller (ref. 10). 
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Features Characterizing Method (Concluded) 


Extending h c and h p (and hence p) to very large values, 
GPC solution approaches LQR solution; hence GPC 
approximates an optimal controller for large p 

Matrices Q and R are control effort and tracking error 
weighting matrices 

Matrix "T (and hence Y ) is invariant if system does not 
change because it is formed from SMP (which are unique) 
Solution for Y obtained by an efficient computational 
procedure 


By extending the prediction and control horizons (and hence p, the order of 
the ARX model) to very large values, the GPC solution approaches that of the 
linear quadratic regulator (LQR). Thus, GPC approximates an optimal 
controller for large p (ref. 5). 

Weighting matrices Q and R are used to limit the control effort and to 
weight the relative importance of the differences between the target and 
predicted responses, respectively. As mentioned earlier, Q and R are usually 
assumed to be diagonal and for Q to have the same value w c along its diagonal 
and R to have the same value w r along its diagonal. The GPC solution offers no 
guarantee of stability and the control weight must be tuned to produce an 
acceptable solution without going unstable. Reducing w c increases control 
authority and performance but eventually drives the system unstable. 

If only the disturbances acting on the system change, there is no need to 
recalculate T (and hence y ) because it is formed solely from the SMP, which 
are unique for a given system (ref. 13). 

The solution for Y described earlier involves forming the matrix products 
yV T and W T . Here, these products are obtained using the computationally 
efficient procedure described in reference 8. 
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GPC Computational Considerations 


System ID in presence of external disturbances for implicit 

inclusion of disturbance model 

Pseudo-inverse based on SVD techniques 

Basis for recalculation of matrices 

On-line/real-time versus on-line/near-real-time 

Batch versus recursive calculation of OMP 

Efficient computational algorithms and coding 


Several considerations dealing with computations should be kept in mind 
when developing algorithms for GPC applications. 

SID should be done with the external disturbances acting on the system 
so that information about the disturbances is embedded in the OMP. However, 
depending on the nature of the external disturbances, it may be possible to 
perform a SID on an undisturbed system and still determine a control law that 
results in satisfactory closed-loop performance. 

The computation of pseudo-inverses should be performed using Singular 
Value Decomposition (SVD) because of the latter’s ability to deal with matrices 
that are numerically ill-conditioned. The use of pseudo-inverses (via SVD) is 
recommended even in cases where the ordinary inverse may seem appropriate 
(such as in the operation (VV ) indicated in the expression for the OMP). 

The basis for recalculation of matrices needs to be selected. For example, 
should the system be re-identified and the control law matrices recalculated in 
every sampling period, every specified number of time steps (block updating), 
or only if some event requiring a re-ID occurs? Should all calculations be done 
in real time, or can some be done in near real time? Should updating of the 
OMP be done in batch mode or recursively? 

Microprocessor speeds are such that it is now often possible to complete 
the full cycle of GPC computations and apply the commands to the actuators 
within one sampling period. However, the need for efficient computational 
algorithms and attendant coding is not expected to diminish. 
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GPC Implementation Considerations 


# Low-pass filtering of measured time histories 

# Cut-off frequency of filter set to Nyquist frequency f N 

® Nyquist frequency chosen so that maximum frequency of 
interest is about 75% of f N 

# Sampling frequency selected to be between 2-3 times f N 

# Length of data stream for SID determined by need for 5-10 
cycles of lowest frequency mode in measured responses 

m Scaling (normalization ) of digitized input and output data 

# Distribution of computing tasks among computers/CPUs 
® Extent of user interaction/involvement 

# Constraints on acceptable values of input and output 


Several considerations must be taken into account when actually 
implementing GPC algorithms in hardware for active controls work. 

The measured response time histories must be passed through a low-pass 
filter with a cut-off frequency f c set equal to the Nyquist frequency /n. The latter 
is chosen so that the maximum frequency of interest is about 75% of f N . The 
sampling frequency fs should be at least twice //v to prevent aliasing. However, if 
fs is made too large the low frequency modes will be poorly identified due to a 
loss of frequency resolution. A sampling rate between 2 to 3 times is 

generally sufficient. Once the sampling frequency has been selected, the 
minimum number of data points that should be used for SID follows from the 
requirement of having 5-10 cycles of the lowest frequency mode in the measured 
response time histories. 

Normalization of the input and output data that is used for SID on the 
maximum actual or expected values of the data is often helpful numerically. The 
procedure will depend on whether the computations are being done in batch mode 
or recursively. 

The computing tasks can be distributed among computers or different CPUs 
on a single computer. The choice will influence the extent of user involvement. 

The values of the input and output data that are being used during closed- 
loop operations must be carefully monitored to ensure that they fall within 
acceptable bounds. This is easily done using IF/THEN-type checks in the code. 
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Illustrative Numerical and Experimental Results 

# Simulations using lumped-mass-spring-dashpot systems 

# Experimental dynamic studies using Cobra stick model 

# Simulated tracking control of XV- 15 in roll maneuver 

# Hover test of stiff-inplane rotor on WRATS 

# Ground resonance test of soft-inplane rotor on WRATS 


A variety of numerical and experimental studies comprised the 
exploratory evaluation phase of the GPC investigation during the period 1997- 
1999. These studies included: (1) Numerical simulations using 3- and 9-DOF 
mass-spring-dashpot systems varying the number and location of external 
disturbances, the number and location of control inputs, and the number, 
location, and type of response measurements (displacement, acceleration); (2) 
Extensive experimental simulations using the BHTI Cobra stick model in which 
the number and location of the shakers used for control and disturbance were 
varied; (3) Preliminary assessment of GPC for simulated tracking response using 
an XV- 15 math model determined using BHTI’s COPTER (Comprehensive 
Program for Theoretical Evaluation of Rotorcraft) flight simulation program 
(ref. 14); (4) A hover test (April 1998) of a stiff-inplane rotor on the WRATS 
testbed to explore the vibration reduction capability of GPC (excitation caused 
by rotor downwash and rpm near resonance); and (5) Ground resonance test 
(Oct. 1999) of a soft-inplane rotor on the WRATS testbed. 
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The 3-DOF mass-spring-dashpot system used in the 3-DOF simulations is 
depicted at the top of the slide. 

The masses were equal to 1.0 and the spring rates were all equal to 1000 
lb/in. The damping coefficients of dashpots 1 and 2 were equal to 0.1 lb/in/sec, 
but the damping of dashpot 3 was set to - 0.1 lb/in/sec to make the system 
unstable (one pair of complex conjugate eigenvalues had a positive real part). 

Control was imposed at masses 1 and 3 (r = 2). Accelerations were 
measured at masses 1 and 2 (m = 2). Sinusoidal disturbances were applied at 
masses 1 and 2 at a frequency f = 6.276 Hz, which coincided with a natural 
frequency of the system. The expressions used to represent the disturbances in 
the simulation are udl = 2*cos(2*pi*f*t) and ud2 = l*sin(2*pi*f*t) and are 
shown in the plot at the lower left. Also, At = 0.05 sec, / = 300, and p = 5. 

Tire oscillatory diverging behavior of the open-loop system is shown in the 
plot at the upper right, which shows the time history of the acceleration 
responses of masses 1 and 2. The loop was closed after 100 time steps. The 
time histories of the closed-loop responses and control commands are shown in 
the plots in the lower right and upper left, respectively. 
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Cobra Stick Model 



In parallel with the extensive numerical simulations that were being 
conducted using the 3- and 9-DOF mass-spring-dashpot systems, bench tests 
were conducted by Bell Helicopter using the “Cobra stick model” shown in the 
slide. This model is a 36-inch long, 50-lb, multiple degree-of-freedom lumped- 
mass dynamic model that approximates the dynamics of a Cobra helicopter. 
The dominant vertical bending modes of the model have natural frequencies of 
about 8 and 23 Hz. Electromagnetic shakers were used in various combinations 
for imposing external periodic and random disturbances, random excitations 
needed for system identification, and the control inputs called for by the GPC 
algorithm. 
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An illustrative experimental result obtained using the Cobra stick model is 
shown in the slide above. The case shown is for sinusoidal excitation of 23.5 Hz 
(near a model natural frequency) at the nose of the model; random excitation at 
the main rotor hub for SID; vertical response at the tail rotor location; control 
input at the hub; f s = 512 Hz; / = 300; m = 1; r = 1; p = 8; w c = w r = 1.0. For the 
case shown here, 300 time points were used to ID the system and imposition of 
closed-loop control was then (arbitrarily) delayed for 50 time steps. 

The top plot shows the time history of the tail boom vertical acceleration. 
The bottom plot shows the time history of the control input during SID, the 
delay, and after closing the loop. The response is reduced dramatically in about 
0.2 sec. 
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Simulated Tracking Control of XV- 15 Tiltrotor 


* Aircraft in airplane mode of flight; 
150 kts airspeed; sea level 

s Control inputs: Collective, pedal, 
longitudinal & lateral cyclic pitch, 

* Responses: Pitch, roll, sideslip, 
and forward velocity 


® COPTER used to generate linearized mass, damping, and stiffness 
matrices for trimmed flight condition and components for B matrix 
® MATLAB used to form A, B, C, and D state matrices and generate 
simulated input-output time history needed for SID/GPC 
® Roll maneuver prescribed 



A cursory examination of using the GPC algorithm as a means for 
actuating the controls necessary to execute a prescribed maneuver was made. 
The subject aircraft was the XV-15 cruising at 150 kts in the airplane mode of 
flight. The controls considered were collective pitch, longitudinal and lateral 
cyclic pitch, and pedal. Responses of interest were aircraft pitch, roll, sideslip, 
and forward velocity. Bell Helicopter’s COPTER flight simulation program 
(ref. 14) was used to compute the linearized mass, damping, and stiffness 
matrices for a trimmed fight condition as well as components for the B matrix. 
MATLAB was used to form the A, B, C, and D state matrices and generate the 
simulated input-output time histories needed for SID and GPC. A one-cycle 
saw-tooth-type roll maneuver was prescribed to be the target response, yr- 
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One cycle of a saw-tooth-type roll maneuver was prescribed to be the 
target response, yr. The time histories of this desired response and the actual 
roll response produced by GPC are shown for comparison at the bottom of the 
slide. The remaining three responses (pitch, sideslip, and forward velocity) had 
negligible magnitudes and are not plotted. 

Tire tune histories of the rudder and lateral cyclic control inputs called for 
by the GPC algorithm to perform the prescribed roll maneuver are shown at the 
top of the slide. The remaining control inputs (collective pitch and longitudinal 
cyclic pitch) had negligible magnitudes and are not plotted. 

Tire primary structural natural frequencies of the aircraft lie between 0.8 
and 10 Hz. The values of pertinent GPC parameters are: w c = 1.0, p = 8, h p = h c 
= 2 p, l = 400, At = .025 sec ,f s = 40 Hz. 
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WRATS Model in Aeroelasticity Branch 
Rotorcraft Hover Test Facility 



The initial evaluation of GPC on the WRATS testbed was conducted in 
April 1998 during a one-week hover dynamics test conducted in the AB 
rotorcraft hover test facility located in a building adjacent to the TDT. A stiff- 
inplane gimballed rotor was employed in this investigation. The model is shown 
above with its rotor blades removed and replaced by equivalent lumped weights 
for a ground vibration test that was conducted prior to the test. Emphasis in this 
initial evaluation of GPC was on active control of vibration using only the 
collective control. To provide a rigorous test of the GPC algorithm, the open- 
loop response of the model was exaggerated by miming the rotor at an rpm that 
nearly coincided with the natural frequency of the wing vertical bending mode. 
Additional excitation of the wing was provided by the downwash associated 
with miming the rotor at a high thrust level. 
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Implementation of GPC Algorithm for 
WRATS Hover Dynamics Test 

Dual-computer approach with user interaction: 



For the hover dynamics investigation, the GPC computations were 
distributed between two computers and user interaction was required to transfer 
the data from one computer to the other. On user command, data required for 
system identification was collected on computer #1 and sent to computer #2 
where SID was performed and the control gain matrices d and ff were 
computed. On user command the control parameters were sent to computer #1 
which used the p latest data sets to (continuously) compute the control 
commands that were sent to the swashplate actuators. If re-identification of the 
system was required, the process was repeated on user command. 

All computations were done using MATLAB on PCs with 500 MHz 
CPUs. 
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Some experimental results obtained during the WRATS hover dynamics 
test that illustrate the effectiveness of GPC in reducing rotor-induced vibrations 
are shown in this slide, which shows the measured time histories of the open- and 
closed-loop responses of the vertical bending moment rear the root of the wing. It 
is seen that the response is dramatically reduced within a second after the control 
is turned on. 
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Implementation of GPC Algorithm for 
WRATS Ground Resonance Test 

Single-computer/dual-processor approach: 


CPU#1 


CPU #2 



For the ground resonance investigation, the GPC computations were 
distributed between two CPUs on a single computer with no user interaction 
required for data transfer as in the earlier hover test. On user command, data 
required for system identification was collected on CPU #1 and sent to CPU #2 
where SID was performed and the control gain matrices d: and // were 
computed. The control parameters were automatically sent to CPU #1 which 
used the p latest data sets to (continuously) compute the control commands to be 
sent to the swashplate actuators. If re-identification was required, the process 
was repeated on user command. 

All the algorithms were implemented in dSPACE software on a PC with 
500 MHz CPUs. 
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Ground Resonance Test of Soft-Inplane Rotor 
on WRATS Testbed 


(8 degrees collective at 75% r/R) 



hpm 


A brief investigation of the use of GPC to actively control the ground 
resonance behavior of a soft inplane tiltrotor was conducted in the AB hover test 
facility in October 1999. For this test, the model blades were modified by 
replacing the stiff inplane flexure at the root of each blade with a spindle 
incorporating a lag hinge and an adjustable viscous damper and spring. The 
open-loop behavior (frequency and damping versus rotor speed) was compared 
with its closed-loop behavior for several values of collective pitch. A GPC- 
based algorithm was used to actively control the cyclic inputs into the (fixed- 
system) swashplate in a manner which produced a whirl of the rotor tip-path- 
plane in the direction and frequency needed to stabilize the critical body mode. 
For the open-loop configurations of the model in which a definitive ground 
resonance instability was observed, use of GPC was found to be strongly 
stabilizing. In particular, damping levels of about 2% critical were noted in the 
rpm range where the open-loop system was unstable. The results for the 
configuration with eight-degrees of collective pitch are shown at the right in the 
slide above. A photo of the model obtained from the video recording made 
during the test is at the left. A pneumatically-actuated snubber system, 
consisting of four horizontal cables attached to the upper end of the pylon in an 
X-pattern, was used to arrest instability and is partially visible in the photo. 
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Concluding Remarks 


Exploratory studies of GPC for active aeroelastic control 
of tiltrotor aircraft completed 

Studies included numerical simulations and experimental 
investigations 

GPC-based MIMO active control system was highly 
effective in increasing stability and reducing response 
GPC appears to be a viable approach for active stability 
augmentation and vibration control of tiltrotor aircraft 

Wind-tunnel tests of WRATS model with GPC needed to 
more fully evaluate method 


Exploratory numerical and experimental studies into the use of GPC for 
active aeroelastic control of tiltrotor aircraft have been completed. A GPC- 
based MIMO active control system was demonstrated to be be highly effective 
in increasing ground resonance stability and reducing vibratory response. While 
these results are quite encouraging with respect to establishing the viability of 
the method, it is recognized that a broader evaluation of the methodology is 
needed to validate GPC-based algorithms for active stability augmentation and 
vibration control of tiltrotor aircraft. 
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Status and Plans 


Initial wind-tunnel evaluation of GPC on WRATS for 
airplane-mode stability augmentation completed 
Method highly effective in increasing stability (damping) 
of critical wing mode 

Continue development and evaluation of GPC-based 
method for active aeroelastic control 

Improve suite of computational algorithms 
Establish robustness to system nonlinearities 
Conduct additional ground and wind-tunnel tests 


The initial wind-tunnel evaluation of GPC on the WRATS testbed for 
airplane mode stability augmentation has recently (April 2000) been completed. 
The method was highly effective in increasing stability (damping) of the critical 
wing mode for all of the model conditions tested. 

Plans are to continue development and evaluation of GPC for active 
aeroelastic control. Emphasis will be on improving the suite of computational 
algorithms comprising the current GPC software system developed for WRATS. 
hi this regard, work is underway on providing for the calculation of the closed- 
loop eigenvalues using closed-loop input-output data and for the recursive 
calculation of the OMP. 

It is expected that work will also be initiated to establish the robustness of 
GPC-based methods to system nonlinearities. Both numerical simulations and 
experimental studies are anticipated. 

Additional ground and wind tunnel tests are to be conducted as necessary 
to evaluate the GPC-based methodology over a broader range of simulated flight 
and operating conditions. 
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